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Abstract 

We consider a three dimensional system consisting of a large num- 
ber of small spherical particles, distributed in a range of sizes and 
heights (with uniform distribution in the horizontal direction). Par- 
ticles move vertically at a size-dependent terminal velocity. They are 
either allowed to merge whenever they cross or there is a size ratio 
criterion enforced to account for collision efficiency. Such a system 
may be described, in mean field approximation, by the Smoluchowski 
kinetic equation with a differential sedimentation kernel. We obtain 
self-similar steady-state and time-dependent solutions to the kinetic 
equation, using methods borrowed from weak turbulence theory. An- 
alytical results are compared with direct numerical simulations (DNS) 
of moving and merging particles, and a good agreement is found. 



1 Introduction 

We consider spherical particles in a viscous flow. The particles move verti- 
cally with their terminal velocity arising from the balance of the gravitational 
effect (fall or buoyancy) and viscous drag. Since, in general, particles of dif- 
ferent sizes rise or fall with different velocities, their trajectories can cross and 
merging can happen. Realistic models of particle merging are quite involved 
and in the present text we are going to consider only two very simplified 
models: either any two particles whose trajectories cross merge, which we 
shall refer to as "free merging" , or merging is restricted to particles of similar 
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sizes (i.e. small particles avoid big ones due to moving along flow streamlines 
bending around the big particle), which we call "forced locality" (deflned in 
Sect. El . 

It will turn out that, although our problem is very simple to state, it is 
very rich in features. The simplifled model can be realized by considering a 
sedimenting kernel in the Smoluchowski coagulation equation. We will derive 
solutions to this equation analytically, and we examine the validity of such 
solutions with direct numerical simulations (DNS), in which we let particles 
evolve individually according to certain rules for collisions and we study their 
overall size distribution. We shall study different stationary regimes, either in 
time t or in the vertical coordinate z, and we will discuss self-similar solutions 
and study the role of local and non-local merging. Whereas time dependent 
solutions of the sedimenting kernel have received a lot of attention in the 
literature [H El [3], the study of height dependence - also treated here - is 
more rare. 

The process we discuss is usually referred to as differential sedimentation 
and has been linked to experimental results and is used to predict rain 
initiation time [HIE]. In particular, the model admits a power law distribution 
consistent with experimental data for aerosols [5] . In our discussion, we will 
obtain this power law as an exact result, rather than by dimensional analysis 
used in previous discussions [U, [7] . We recognize this result as a Kolmogorov- 
Zakharov (KZ) cascade of the volume integral, similar to the solutions that 
arise in wave turbulence. Solutions to the coagulation equation with a KZ 
cascade have been studied in general [H Ej, and with a kernel describing 
galaxy mergers in particular [TO]. 

We flnd that in the free-merging model the locality assumption necessary 
in dimensional analysis and the KZ spectrum fail to hold We will obtain 
an analytical solution for such a non-local system, and verify this with DNS. 
We will study self-similarity for both the forced-locality model and the free- 
merging model. We will perform DNS for inhomogeneous solutions that are 
self-similar in the spatial variable z. 

The starting point of our analysis is to write a kinetic equation for the 
coagulation process in Sect. 12. 1[ In Sect. Owe flnd the Kolmogorov-Zakharov 
solution for the kinetic equation. Sect. H] discusses the dominance of non-local 
interactions in the system. We study self-similarity of our model in Sect. [5], 
and we analyze locality of such solutions in Sect. El where we present numer- 
ical data. Finally, we introduce a "super- local" model in Sect. [TJ reducible 
to Burgers equation. 
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2 The model 

Let us denote by a the volume of a spherical particle and by r its radius, 

a = nr^ , k = Air/S . (1) 



u(r)=(2gp/9T]) r 




Figure 1: A particle's terminal velocity u is determined by its radius r. 
Larger particles will have a larger terminal velocity, depicted by the arrows, 
following definition ([2]). (Created by T.H.M. Stein) 



The Stokes terminal velocity of a rigid sphere of radius r with no slip 
boundary conditions is given by the formula fHX 0, [T2] 



cr 



"^gjPf - Pp) 
9rif 



(2) 



where g is the free fall acceleration, pf and Pp are the density of the sur- 
rounding fluid and the particle respectively, and rjf is the dynamic viscosity 
of the surrounding fluid. 

Experimentally, the formulae ([2]) are valid for air bubbles in water at 20°C 
with r < 1mm, and these bubbles can be considered spherical. Slip-flow 
corrections can be necessary for other gases and fluids [12]. The following 
data for water droplets and particles in the atmosphere can be found in 
Pruppacher and Klett |5J. For droplets, corrections to ([2]) are necessary 
when r > 30yum, which changes the formula's dependence on r^. They can 
be considered spherical for radii up to 535/im. For atmospheric particles, ([2]) 
can be considered to depend on for large particles. However, atmospheric 
particles are generally not spherical and will thus require other corrections. 
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Despite physical complications, we will assume and ([T]), and we will 
express both in terms of volume a, 

r{a) = /c-V3^i/3 ^ ^(^) = c«:-2/V2/3 . (3) 

We compute this model using direct numerical simulations in a periodic 
box of 10 X 10 X 10cm with particles that are defined by their x-, y-, and 
^-coordinates and by their volume a. At each time step the particles move 
according to their fixed terminal velocity, using definition ([2]). We fix our 
parameter c such that a particle of radius 0.1cm moves upwards with velocity 
20cms~^, which resembles the situation of air bubbles in water [12]. 

The particles are generated at a range of small cr, with their smallest 
volume (Jo ~ 4.2 ■ 10~^cm^, equivalent to a radius r = 0.01cm. They are 
removed from the system once they become larger than IO^ctq, or r ~ 1mm 
and are assumed to be spherical at all sizes for computational purposes. 
With different velocities, the particle trajectories may cross, and depending 
on the rules of interaction they can then merge. These rules are governed by 
collision efficiency, which will be explained in Sect. 12. 1[ 

2.1 The kinetic equation 

We suppose that the distribution of particles can be adequately characterized 
by density n{a, z, t) (the number of particles of volume between a and a + 
da, per fluid volume V per da, at the vertical coordinate z and at instant t). 
In particular we suppose here that the dependence of particle distribution on 
the horizontal coordinates can be averaged out. This hypothesis is valid if the 
dynamics do not lead to strongly intermittent distribution in the horizontal 
directions, for example if the fluid is well mixed in the horizontal directions. 
Our numerical simulations appear to support such a mean field approach 
well, and in future work it would be interesting to examine theoretically why 
this is the case. 

The goal of this section is to derive a kinetic equation for n - also called 
Smoluchowski coagulation equation [13] - using a kernel describing differen- 
tial sedimentation. We write the collision integral, which expresses simply 
the fact that two particles of volumes ai and (J2, with ai + a2 = c, can merge 
to give a particle of volume a (inflow), or a particle with volume a can merge 
with any other particle of volume ai > and give a particle with volume 
(72 = cr + o"! (outflow). Also, wc determine the cross-section of interaction 
between two particles by the condition that particles merge upon touching, 
that is if their centers are at a distance at most ri + r2, which gives the geo- 
metric cross-section of 7r(ri + r2)^. Finally the collision rate between particles 
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of volume cxi and (T2 is taken to be proportional to their relative velocities 
\u{ai) — M(cr2)| and to their number densities rii and n2, which is a mean field 
type hypothesis. 

The left hand side of the kinetic equation contains the advection term 
dtn + udzU, which we shall also denote as the total derivative dn/dt, while 
on the right hand side we put the collision integral. Note also the shorthand 
n = n{a,z,t), u = u{a), rii = n{ai,z,t), Ui = u{ai), ri = r{ai) and similar 
for n2, U2 and r2. Thus we find 

dtTi + udzU = (4) 
1 

+ - dcridcr2 |m2 - ■Ui|7r(ri + r2)^nin2(5((T - (Ti - 0-2) 
^ Jo 



daida2\u — n2|7r(r + r2Ynn25{ai — a — (T2) 
do"ido"2|^ — Ui\iJ:{r + rif'nni5{o2 — o — a\) 



It is useful to express the u and r in terms of o using ^ , 

dtn + cKT'^l^o'^l^dzn = (5) 



1 2/3 






1 l/3\2 r/ 
+ (T2 ) 721^20(0- 


- ai - 


- 0-2) 








V3\2 e/ 

+ (T2 ) nn20[ai 

V3\2 


— 0" — 

— (T — 


0-2) 

^1) 



Let us introduce the interaction kernel i^(cri,cr2), 

—4/3 

rv-/ N CK I H 2/3 2/3|/ 1/3 , 1/3n2 fa\ 

K{ai,a2) = |(t/ -a/ I (a/ + ^2 ) , (6) 

which for a general kernel K reduces Eq. (jl]) to the Smoluchowski equation. It 
is useful to note that our kernel (Ej) is homogeneous in a, with K{(ai, (a2) = 
C'^/^i^((Ti, (T2). We also introduce the collision rates 

Rai2 = (y2)nin25{a - ai - 0-2) (7) 

with Ria2, R2ai defined analogously. Now the RHS of Eq. can be written 
in a compact form 



dn 
dt 



p+oo r+00 

/ ddi / d(T2 {Ral2 - Rla2 " -^2,71 

Jo Jo 
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2.2 Characteristic timescales 



We study the physical relevance of Eq. ([5]) by comparing its characteristic 
time Tds with the characteristic residence time in a typical system, Tg = L/u, 
where L is the vertical extent of the system, and u is as in Eq. ([3]). To find 
Tds, we note that n ~ ^ and we introduce the volume fraction f ~ so 
that: 

V 

n — . 

Now, using the kinetic equation ^ we can write 



1 _,1 0^0/Q^0/Q_1 

72 



Thus we find the following relation between the characteristic times: 

= -2/3 2/3 - ' ^0 

where we recall that cr^/^ = /t^^/^r and approximate k^^/^tt 2. From [5J we 
find that for a cumulus cloud, typically L ~ lO^m, r ~ lO^^m, and v ~ 10^^. 
Thus, we find that Tg/rds ~ 10^, which implies that the kinetic equation 
is relevant in a cloud system with gravity when we regard time and length 
scales. 



2.3 Collision efficiency 

The kinetic equation ([5]) allows merging of particles of any sizes, without 
any discrimination. We shall refer to this "free merging". More re- 

alistically one should also take into account the collision efficiency between 
particles. We define collision efficiency £12 = S{(7i,a2) between particles of 
volumes ai and (T2 as a number between and 1, which enters the collision 
integral by multiplication with the collision rates R, so Rau would be re- 
placed by Rau^u and more generally for example the integrand of Eq. ([H]) 

would become RaU^U - Rla2^a2 - R2al£al- 

In particular, one could restrict merging to particles of similar sizes, tak- 
ing into account that small particles cannot collide with much larger ones 
because they bend around them along the fluid streamlines. In the simplest 
such model which will be considered later in this paper, 

^12 = 1^ if 1/g < f^i/f^2 < g, ^^^^^ 

1 otherwise. 
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Figure 2: Without applying the efficiency kernel S, particles merge whenever 
they cross. Including £ with small g, only situation B is allowed, i.e. only 
particles of similar size may merge; particles of dissimilar size (situation A) 
are allowed to cross one another without merging. (Created by T.H.M. Stein) 



where g > 1 is the number representing the maximal volume ratio for the 
particle merging. Compared to a more involved form of collision efficiency 
used by Valioulis et al. the simplified kernel we use mimics the behav- 
ior for particles with r = 0.01cm which is similar to the regime we study 
numerically. We will refer to the model with finite q as "forced locality" . 



2.4 Scaling argument 

For our simple setup one could derive a steady state solution merely by 
physical and dimensional arguments, following Friedlander [IS], Jeffrey [7], 
and Hunt The main remark is that at steady state, the system has a 
constant fiux of volume. The total volume of particles per unit volume of 
fiuid that passes from particles smaller than a to particles greater than a is 
of the order: 

We can estimate from the kinetic equation ([8]) and equations (JTj) and ([6]) that 
dn/dt ~ a^R, with R ~ Kn'^a~^ and K ~ a^^^. If we assume that n ~ a'^, 
we find that dn/dt ~ a^/^^^'^, and we obtain the scaling cr^^/^+^^ for the 
volume fiux (fT^ . For constant fiux, we arrive at z/ = —13/6, or n ~ a^^^^^. 
Naturally, the dimensional analysis assumes locality of interactions. 
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3 Kolmogorov-Zakharov solution 



One of the simplest questions one can ask with respect to the kinetic equa- 
tion ([5]) is if it allows for a scaling stationary solution of non-zero flux. 
Such a solution, if one exists, is called a Kolmogorov-Zakharov (KZ) spec- 
trum because, like in the classical Kolmogorov spectrum, it corresponds to 
a cascade of a conserved quantity (total volume occupied by particles in our 
case) [Hllin]. In this section we investigate the scaling exponent and existence 
of such solutions. 



3.1 Zakharov transform 

A derivation of the KZ solution can be achieved through the technique of 
the Zakharov transform Ti5J. Let us consider a steady state (i.e. time and 
space independent) solution of Eq. ([5]) of form n ~ a" , and let us aim to find 
V. Note that this is a reasonable thing to look for, since we can easily see 
from Eq. ^ that our collision integral is a homogeneous function in the a 
and in the n. 

We start by expanding our collision rates from equation (JTj) using equation 
(EI), and obtain the following equation in a: 

—4/3 

' Vr 2/3 2/3|/ 1/3 , 1/3n2 v i/r/ \ 
Ral2 = ^ ^2 -^1 1(^1 +(^2 ) fTia2(5(a-ai-(T2) 

where Ria2 and R2ai are expanded similarly. We then continue by non- 
dimensionalising the rates R by writing ai as cr^o" and (T2 as cTgcr, so 

R.12 = (13) 



—4/3 

'^^ ^ ^1/3+2;^ 1^/2/3 ^'2/3,/ /1/3, / 1/3x2 /V / /■ 
(y ' ^2 -^1 K^l +0^2 ) ^1 ^2 t>(l - ^1 - Cr2, 



and Ria2 and R2ai are transformed in a similar way. 

The Zakharov transform consists in passing in -Rio-2 to new variables 5"i 
and a2 defined by 



, 1 ^' ^2 

CTl CTi 



This way, we obtain 

Rl.2 = (14) 

—4/3 

CK TT 2z/+l/3~-l/3-2i/|~2/3 ~2/3|/~l/3 , ~ 1/3x2 ~ i/ ~ i/ r/ . ~ ~ \ 
^ <^ ' ^2 -^1 1(^1 ^ o4 ) a^O^b^l- Ox- 02) ■ 

A similar expression is derived for R2a\- 
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Combining the transformed terms and dropping primes and tildes, we 
transform the compact kinetic equation ([H]) 

0= da J da, (1 - a-'"'-'^ - a^,'"''"')R.,, . 
Jo Jo 

Here, we note that the integration variables for Ria2 become daidcr, = 
cr^ajf^dcTido",, with a similar transformation in R2ai- Now, if we choose u such 
that —10/3 — 2z/ = 1, then we have the factor 6{1 — o"i — o'2){l — (Xi — cr,) = 
appearing in the integrand, which solves the equation, i.e. u = —13/6 is the 
candidate for the KZ exponent. This method of derivation can be applied to 
various kernels for the Smoluchowski equation 

Let us note that our exponent u is that of n{a). In literature, one com- 
monly finds the radius distributions, n(r), which can be expressed in terms 
of n{a) from the relationship n{a)da = n{r)dr. Thus, n(r) = n{a)da/dr oc 
^3u^2 ^ ^3u+2^ g^^^ therefore = 3u + 2 = -9/2 |7j. 

However, the KZ spectrum is only a true solution of Eq. ([5]) if the collision 
integral on the RHS of this equation (prior to the Zakharov transformation) 
converges. This property is called locality, and it physically means that the 
particle kinetics are dominated by mergings of particles with comparable 
(rather than very different) sizes. Convergence of the collision integral on 
general power-law distributions will be studied in Appendix \M We will see 
that (without modifying the model to enforce locality) the —13/6 scaling 
exponent gives rise to non-local interaction between the particles both with 
the smallest and the largest particles and, therefore, the KZ spectrum is not 
a valid solution in this case. 

3.2 KZ spectrum in the system with forced locaUty 

Locality of interactions, and therefore validity of the KZ solution, are imme- 
diately restored if one modifies the model by introducing the local collision 
efficiency kernel as in definition ffTTj) . This kernel is a homogeneous func- 
tion of degree zero in a and, therefore, the KZ exponent obtained via the 
Zakharov transformation remains the same. In Fig. [3] we can see that the 
Kolmogorov- Zakharov scaling appears in a system with forced locality. 

4 Kinetics dominated by non-local interac- 
tions 

As an alternative, we may assume that the dominant interactions are non- 
local and find a cut-off dependent stationary solution. This is relevant if it is 
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Figure 3: Distribution of particle volumes averaged over several times after 
140,000 time steps for the forced locality situation with q = 2. The dashed 
slope represents the —13/6 KZ spectrum (compare with [11]). 



not desirable to use the collision efficiency models which guarantee locality 
(for instance using the kernel ffTTl) ). In this case one should accept the fact 
the kinetics are dominated by non-local interactions, and that the low-cx 
or/and high-cx cut-offs dominate the collision integral. In fact, such a non- 
locality can allow us to significantly simplify the kinetic equation and reduce 
it to a differential equation form. As shown in Appendix El contribution to 
the collision integral from non-local interactions with the smallest particles 
((Ji <^ a) is 



where we have dropped the explicit dependence of the upper integration 
limit on a, since the integral is divergent as (Jmin — ^ (this is the hypothesis 
of non-locality), so the dependence on the upper bound is a sub-dominant 
contribution. 

The contribution to the collision integral from non-local interactions with 
the largest particles (ai ^ a) is 





(15) 



where 




(16) 
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Similarly to above, here the lower integration bound is omitted. 

Putting these two formulae together, we obtain the following effective 
kinetic equation for the cases when the non-local interactions are dominant, 

^ = -cMa'/'n) - c^n , (17) 

where constants Ci, C2 are defined in the formulae (fT5|) and (fT6|) . Note that 
this equation fll7l) is valid when the non-local interactions with the smallest 
and with the largest particles give similar contributions, as well as in cases 
when one type of non-locality is dominant over the other. 

In steady state dn/dt = and the solution of the resulting ordinary 
differential equation is 

n = Co-'l^'e^'''"" , (18) 

with C being an arbitrary positive constant. Note that the constants C and 
C2/C1 appearing in the solution ffTSl) can be related to the "physical" data of 
Cmin, cTmax and n((Tmin), through Eqs. ( |T5l) . ( |T6|) and ( |T8|) . We obtain 



exp 

n{a) = n(o-min) 



-1/3 

log^ 



4/3 



(19) 



The solution ffT^ is interesting since it is not a pure power law. For 
large a we have n ~ Ca~^^^ which is a limit when absorption of the smallest 
particles is much more important than being absorbed by the large particles, 
i.e. when the first term on the LHS of Eq. ( fTSi) is much greater than the 
second one. This limit corresponds to a cascade of the number of particles 
(not their volume!) which is a conserved quantity in this regime. 

In Fig. m we show our numerical results for the non-local model. Particles 
are produced uniformly in space with volumes ranging from (Tq to Suq, and 
particle density within this size range is kept constant in time. Particles 
are removed from the system once they reach cXmax = 10^o"o, with probability 
p{a) = 1— exp~"('^^'^""="') with a <^ 1. The original results have been averaged 
over neighbouring data points to obtain the continuous graph in Fig. |H We 
also used Eq. ( fT9l) and find that with appropriate parameters this solution 
fits the numerical data. 

We can check our hypothesis of dominance of non-local interactionsdi- 
rectly by counting the number of collisions within a certain timeframe at 
statistical steady state. Namely, for each size bin we count the number of 
collisions leading to a particle entering the bin, and the number of colli- 
sions leading to a particle leaving the bin. We distinguish between local and 



11 




Figure 4: Averaged distribution of particle sizes for the situation without 
forced locahty ("g = oo") after 200,000 time steps. The vertical dotted lines 
bound the inertial range at cTjam = 3cro and o"max = IO^ctq- The dashed curve 
represents the fit conform eq. ffT^ . with a^i^ and a^aax given by the bounds 
of the inertial range, and n(amin) = 1.5 • 10^°; the dash-dot slope represents 
a power law of a^^^^. 



non-local collisions using the particle size ratio q*, i.e. if 1/10 < q* < 10 
we consider the collision local, and non-local otherwise. For non-local colli- 
sions, we distinguish between a collision with a very large particle and a very 
small particle. In the kinetic equation ([5]) (which we do not rely on in our 
procedure) this would correspond to splitting the collision integral as follows: 

— =+ da-i/((7i,cr - CTi) - / d(Ti/(o-i,cr) 

+ / d(Ti/((7i,(7-(7i) - / dai/(ai,(T) (20) 

"^max 

dcTi/ (0-1,0-) 

go- 



where 



/(cri,o-2) = K(o-i, 0-2)721712 . 

We perform DNS and for each collision that occurs we count its contri- 
bution to the different collision regimes as mentioned above. Our results 
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10° lo' 10^ ct/cto 



Figure 5: Number of collisions N per bin [LI'^ctq, 1-1'^'''^o'o] over 10,000 time 
steps, which lead to a particle entering or leaving the bin. Triangles: con- 
tribution due to collisions with large particles; circles: contribution due to 
collisions with similar sized particles; squares: contribution due to collisions 
with small particles. Filled and open symbols correspond to number of par- 
ticles entering and leaving the bin respectively. 



are shown in Fig. [5l We notice that once collisions with small particles are 
counted at cr/ao = q, with g = 10 in this figure, their contribution dominates 
almost immediately, and remains dominant for the entire inertial domain. 
We can also see that collisions with larger particles are only dominant in 
the forcing range a < Suq, and collisions with similar sized particles only 
marginally dominates in the intermediate regime for Sctq < a < SOctq. 

5 Self-similar solutions 

KZ solutions studied in Sect. [3] are valid stationary solutions of the kinetic 
equation ([5]) in the systems modified by introduction of a local collision effi- 
ciency (e.g. using the model ( fTTl) ). We have argued in Sect. Hlthat without 
such an enforced locality the non-local interactions are dominant which re- 
sults in a prediction for the steady state given in Eq. ffT^ and which is 
qualitatively confirmed in direct numerical simulations of the dynamics of 
particles. 
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However, both of these approaches assume homogeneity in space as well 
as a sink at large volumes (i.e. removing particles from the system when they 
reach a certain large size). These two conditions cannot be made realistically 
consistent because there is not a physical mechanism that could remove large 
particles from the bulk of the fluid. 

Thus, it is more realistic to consider one of the following solutions: 

• time-dependent, height-independent solutions without a sink 

• height- dependent, time-independent solutions with a sink at a given 
height (i.e. for bubbles in water an interface with air at a given maxi- 
mum value of z). 

Both situations can be described by self-similar solutions of the kinetic equa- 
tion In the following derivations of the self-similar solutions we will 
suppose locality, in the sense that the dimensional analysis leading to the 
results supposes no dependence on the cut-off scales (Tmin and (Jmax- Validity 
of the locality hypothesis will have to be examined a posteriori. 

We will start by considering the particle model without forced locality, 
and later we will proceed by adding the effect of local collision efficiency 
followed by a super-local model leading to Burgers equation. 



5.1 Height dependent solutions 

Let us start with the analysis of the time-independent state. We look for a 
solution n that is self-similar in the sense that it verifies the scaling relation 

n{a,z) = z^'hiz'^a) . (21) 

To determine the exponents a and f3 we need two relationships. The first 
one is that Eq. should give an equation on h as follows: introduce the 
self-similar variable r = z^a to replace all occurrences of a, then Eq. (jS]) can 
be written as 

T^'^z'^-i^-^ahiT) + ^Th'ir)] = z^^-lP / dn / dr, (T.12 - - T^ri) 

Jo Jo 

(22) 

with the rate 

—4/3 

Tru = '-^\rr - rriirl^' + r'/rHn)HrMr - n - r,) 

with Tit-2 and T2ri defined accordingly. We need to have equal powers of z 
on both sides, which gives 

2 7 
a--p~l = 2a--P . 
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The other relationship expresses constant flux of mass through a given 
height z. Since droplets of volume a move with speed u = u{a), this flux 
is / n{z, a)uada. With h and r this becomes / z°'h{T)z~'^^^^T'^^^z~^Tz~^dT. 
The total power of z should be for z to vanish from this expression, which 
gives us the second relationship 

a-^P = 0. 
Combining the two relations on a and (3 we find 

« = -^, /5 = -l, (23) 

implying 

n{a,z) = z-^/^h{a/z) . (24) 



5.2 Time dependent solutions 

Let us consider a self-similar distribution independent of z but dependent on 
time, of the form n{a,t) = t'^h{t^a), where t = t* — t and t* is a constant, 
the meaning of which will become clear shortly. The left hand side of Eq. ([5]) 
is replaced by dtu = at'^~^h{t^a) + i3t°'~^^~^ah'{t^a). Upon introducing r = 
t^a, this becomes t"~^[a/i(r) + jSrh^r)]. The right hand side of Eq. fl22|) is 
unchanged except for replacing z by t. We thus obtain our first relationship 

1(3 - a = 1. (25) 

One could think that the second relation should come from the conservation 
of mass / n{t, a)ada = J t"h{T)t~^Tt~^dT. However, this condition is in- 
correct because the self-similar solution in this case gets realised only in a 
large-cr tail whereas most of the volume remains in the part which is not self- 
similar. This situation is typical of systems with finite capacity distributions, 
and it has been observed previously for the Alfven wave turbulence ^17j and 
for the Leith model of turbulence [IH]. Thus, we have 

n{a, t) = it* - tyh [ait* - t)3("+i)/7) . 

As in the case of the Alfven wave turbulence [17], it is very tricky to establish 
how to fix the second constant a but it can be found via numerical simulations 
of the kinetic equation 

The above self-similar solution describes creation of infinitely large parti- 
cles in finite time, which rise with infinitely large velocities. Thus, no matter 
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how large our system is, close to the moment t = t* there will be parti- 
cles that travel across the entire height in short time and, therefore, the 
z-independency assumption will fail. Note however that even close to the 
singularity moment t = t* the total volume fraction of such large particles 
remains small. We will study further details of such self-similar solutions 
using the "super-local" model in Sect. 17.21 

6 Locality of the self-similar solutions 

Locality of interactions was assumed in the derivation of the self-similar so- 
lutions in Sect. 15.11 This does not need any further justification if a local 
collision efficiency like in Eq. (fTTj) is used. However, in the case of cut-off 
free interaction kernels that assumption needs to be verified. In order to 
examine its validity we will now establish the asymptotic behavior, at small 
r and at large r, of the self-similarity function h^r) introduced in Sect. [5l 
We shall make the hypotheses (to be verified below) that at very large r the 
collision integral is dominated by contributions of the range of much smaller 
T and, conversely, that at very small r the collision integral is dominated by 
contributions of the range of much larger r. 

Let us start with the large r case. Under the assumption for this range 
that we formulated in the previous paragraph, the distribution in this range 
evolves as in Eq. (|T3|) . i.e. in the 2;-dependent steady state we have 

udzH = —Cidcr{(j'^^^n) , 

which for /i(r) reduces to 

r2/3[a/i + prh'] = -cy^^'lh + rh'] . 

o 

Both sides are homogeneous in r, but the left hand side is of degree 1/3 
higher than the right hand side, so its dominant contribution should cancel, 
leading to the asymptotics /i(r) ~ r~"/^, and substituting values of a and 
P from Sect. 15.11 we get /i(r) ~ r~®/^. According to the results summarised 
in Table [11 such —8/3 tail corresponds on one hand to convergence of the 
collision integral at the large a limit (as assumed in the self-similar solution) 
and, on the other hand, it corresponds to dominance of interactions with 
much smaller r's as was assumed for derivations in this section. 

Let us now consider the small r range. As we have hypothesized above 
about this range, the dominant contribution to the collision integral now 
comes form the non-local interaction term with large particles, which for 
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small a behaves as given in Eq. ffTUj) . leading to 

udzTi = —C2n , 

which for /i(r) reduces to 

T^^^[ah + pTh'] = -C2h . 
This can be solved explicitly and yields 

h{T) = Coe^^-^^V-°/^ = Coe-^-^^V-«/=^ , (26) 

where Co > is an integration constant and the last member has values of 
a and /3 substituted from Sect. 15. 1[ Thanks to the very strong stretched 
exponential decay of h at small r the self-consistency of our hypotheses is 
straightforward to verify. At the same time, such fast decay at small r ensures 
convergence of the collision integral at the a = limit. 

We have therefore proven that our self-similar solutions are local. Note 
that this result is remarkable because, in contrast with the KZ solution, the 
locality property holds even without introducing a local coUisional efficiency 
factor. 

6.1 Numerical verification of the height dependent so- 
lutions 

We have performed direct numerical simulations of the set of particles cor- 
responding to the set-up where one should expect the self-similar behav- 
ior. Namely, we generate particles with distribution n{a) = sin(7r((T — 
a^) /I3)a^'^^^ and with vertical coordinate < 2; < 0.5 and we take them 
out of the system as soon as their center has crossed the surface at z = 10. 

The results for the simulation with free merging are shown in Fig. [6l 
A rescaling to self-similar variables has already been done. We see that 
profiles at different z collapse, which confirms the self-similar character of our 
distribution with the self-similarity coefficients a = —8/3 and (3 = —1 found 
in Sect. 15.11 Moreover, we observe that our profile at large r is consistent 
with the —8/3 power law found above. 

We have also performed computations with the forced locality model as 
given in Eq. (ITT!) with q = 2. It comes to no surprise that the observed 
distribution is also self-similar (since the assumed locality has become even 
stronger). Naturally, the shape of the self-similar function /i(r) is now dif- 
ferent. It is interesting that instead of the —8/3 scaling we now see a —5/3 
slope. We will see in the next section that such a slope can be predicted by a 
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Figure 6: Distribution of particle volumes after 39,000 time steps for the 
situation without forced locality ("g = oo"). The graph is presented in self- 
similar variables according to Eq. flM|) . The markers identify the spectrum 
for z = 1.75 (x); z = 3.75 (o); z = 5.75 (+); z = 7.75 (*); z = 9.75 (0). The 
dotted slope represents a -8/3 power law. 



"super-local" model where the integral kinetic equation is replaced by an 
effective differential equation preserving the scalings of the local interactions. 
In the range of large r we observe an exponential decay /i(r) ~ exp(— 6r) 
(where 6 is a constant), see Fig. [71 As will be shown below, these results are 
also predicted by a (regularised) "super-local" model. 



7 Burgers equation for local interaction case 

We will now study the systems with forced locality in greater detail by in- 
troducing a "super-local" model which preserves the essential scalings of the 
original kinetic equation 1^, i.e. 

dtn + ud,n = -a-^d^{a^^/^n'^) . (27) 

Particularly, Eq. fl^Tj) has the same self-similarity exponents as those found 
in Sect. El in either case of height dependent or time dependent self-similar 
solutions. We see that on the right hand side n appears squared, making the 
equation reminiscent of Burgers equation. We are going to pursue this idea 
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Figure 7: Distribution of particle volumes after 23,000 time steps for the 
forced locality situation with q = 2. The graph is presented in self-similar 
variables according to Eq. [2H The markers identify the spectrum for z = 1.75 
(x); z = 3.75 (o); z = 5.75 (+); z = 7.75 (*); z = 9.75 (0). The dotted slope 
represents a —5/3 power law, and the dashed curve shows Ar~^/^ exp~'>'^, 
made to fit the data at r = 6. 



below, by studying the simpler cases of stationary solutions of this equation, 
either in z or in t. 



7.1 Height dependent solutions 

If we look for steady state in t only, then Eq. ( 1271) reduces to 

We turn this into Burgers equation by introducing new variable s such that 



a = s' 



and the new function 

g{s) = As'-niais)) . 

Then d^g = -{AX)-\s^'-^^/^+^^sis^^^/^-^^'g^). If we set fi-8X/3 + l = and 
13A/3 — 2fi = and (AX) = 2 then we recover Burgers equation: 

dz9 = -gdsg . (28) 
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This happens for A = 2, /i = 13/3 and A = 1. 

Conservation of total particle volume leads to the conservation of the 
integral / g{s)ds, and we deal with the usual Burgers dynamics even for the 
weak solutions (i.e. any regularisation of this equation should conserve the 
volume). In this case we get no finite-time singularity since A and A are 
positive. We will use the analogy of (1^ with Burgers equation and assume 
a discontinuity in our function g would be a shock in the equivalent Burgers 
system. The sawtooth shock can be seen to evolve such that at "time" z the 
shock is at s* ~ z^^'^ and its height is ~ z~^^'^ (hint: write ds^/dz = g^/2 
and = B where B is constant). For the original variables this gives 
cT^, ~ z'^^'^ = z and uq ~ = z~^^^. One then sees that this solution 

is self-similar with the scaling we have found above. In fact 

'2-8/3(a/^)-5/3 ifa<2, 
if a > 2;. 

Remarkably, the —5/3 scaling of the self-similar function h{r) is indeed 
observed in the numerical simulation of the particles with the forced locality 
collision efficiency, see Fig. [71 This fact indicates that, in spite of simplic- 
ity, the super-local model ( 127]) is indeed quite efficient in predicting certain 
essential features of the particle kinetics. However, we have not observed 
any signature of a shock in our numerical results. Such a shock should be 
considered as an artifact of super-locality which is smeared out when a finite 
interaction range is allowed. 

In fact, following the method exposed in Sect. 4.2 of ref. p], it is also 
possible to obtain the asymptotic behaviour of n{a, z) for large t = cr/z 
(see Sect. 15.11) . This is beyond the reach of the Burgers model 0. Follow- 
ing ref. [2\ and using notation from our Sect. 15.11 we introduce the ansatz 
/i(r) ~ AT~^e~'^'^ , where A, 7 and 6 are real constants, of which we shall only 
determine 6 here. With this ansatz and using the flux formulation described 
in Appendix [B| in particular Eqs. (!29l) and (130|) . we can write Eq. (122|) as 
(note that we take the values of a and (3 from Eq. 



T^/^[-lAr-'e-^^ + {9- ^T)AT-'e-^^] = 

T-^dr / dri / dT2 K{n,T2)A^rl-^T2^e~^^^'+^'^ 

Jo Jt-ti 

The left hand side scales as r^/^-^'g-T-r -ftrhile the right hand side can be 
seen to scale, for large r, as r^/3~2^'g-7T order to see this, note that 



^Even if we added diffusive regularization to the Burgers model to account for not 
strict super-locality, we would get the incorrect z^^^^ exp(— 7(t/z) behaviour, where 7 > 
is some constant (see also Appendix [B|) . 
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g-7(n+-r2) attains its maximum over the integration domain along the segment 
Ti +T2 = T with Ti,T2 > and becomes much smaller for T1 + T2 — T > 7"^, so 
that the effective integration domain is a band of width of order 7"^ around 
the segment T1+T2 = t). In order for the two sides to have the same scaling we 
must have 9 = 2/3. Then /i(r) ~ Ar-^/^e"^^ and n(cr, z) ~ Az-^a-^^^e-^"/' . 

7.2 Time dependent solutions 

Let us now seek 2;-independent solutions of Eq. ( |27j) . In this situation the 
latter reduces to 

dtn = -a-'d^ia^^^^'n^) . 

We turn this into Burgers equation as above, introducing s and g{s) as above. 
Then dtg = -{AXyh^'-'^^+^dsis^^^^^-'^^'g'^). If we set /z - 2A + 1 = and 
13A/3 — 2/i = and AX = 2 then we recover Burgers equation. This happens 
for A = —6, fi = —13 and A = —1/3. 

In order to know what happens at shocks we need to know what quan- 
tity is conserved by evolution, even at shocks. We know that the original 
system conserves the volume J nada, which translates for g to conservation 
of (X/A) J g{s)s'^^~'^~^ds, and since 2A — /x — 1 = this simply means con- 
servation of J g{s)ds. Thus once again we really deal with the usual Burgers 
dynamics. 

If the initial distribution of n is peaked around cxo with height no then the 
initial distribution of g is peaked around Sq = with height = AsqUq. 
It is convenient to suppose that the peak is of compact support, say between 
(Ji < (72, corresponding to si > S2- Since n (the particle density) is positive 
but A is negative, g will be negative and shocks will move towards smaller 
s. The peak evolves to give a shock, which will have formed at some s > S2- 
To good approximation we get a single sawtooth shock which moves towards 
and reaches it in finite time, which for n means (since A < 0) that there is 
a finite-time singularity at infinite volume. 

The important feature is that the shock in g will arrive at s = at some 
finite time t*, and for t close to t* its height and speed are approximately 
constant, say height g* and position s = tw* where t = t* —t. This translates 
for n to a jump of height A~^s~'^g* = A~^{tw*)^'^g* oc t^^ at position a = 
s'^ (X t^. This is compatible with self-similarity n{a, t) = f^hifl^a) only for 
exponents a = —fi = 13 and j3 = —X = 6, which satisfy the condition from 
Eq. m- 

Note also that, since g can be considered to be approximately constant 
behind the shock (i.e. towards large s) , the distribution of n behind the 
jump (i.e. towards small a) is like a~^^^^, which is a finite capacity power 
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law, as required by conservation of total initial finite mass. 

Since self-similarity only appears in the tail of the distribution, and the 
tail has finite capacity, it is difficult to obtain good statistics in numerical 
simulations for this model. In the tail, there will be very large particles, but 
the void fraction will be large too, as J nada is constant, resulting in a sparse 
data set in the numerical simulation. 

8 Concluding remarks 

As we have seen, the very simple model in which particles move at their 
terminal velocity and merge upon collision appears to be very rich in features. 
For this model, we have derived the Smoluchowski kinetic equation with 
a kernel for differential sedimentation. 

First of all, we considered a setup analogous to one used in turbulence 
theory where small particles are produced and large particles are removed 
from the system with a wide inertial interval in between these source and sink 
scales. We obtained a KZ spectrum (Fig. [3]) and showed that it is relevant 
for the systems with forced locality but irrelevant in the free-merging case. 
In the latter case we derived a model ( ITTj) in which the dominant interactions 
are non-local and we obtained its steady state solution in Eq. ( |T8l) . which 
was verified with DNS (Fig. H]). 

We have also considered self-similar solutions which are either height 
dependent or time dependent. This was done for both the kinetic equation ([5]) 
and for a model with "super-local" interactions (1271) . For the time dependent 
dynamics, we predicted a finite-time creation of infinitely large particles. The 
solutions for height dependent dynamics were verified with DNS. Although 
most particle distributions in the atmosphere are height dependent [3], the 
relevance of self- similarity in such distributions requires further study. 

Our theoretical results were obtained from the kinetic equation ([5]) which 
is essentially a mean field approach. Thus, it is intriguing that such theo- 
retical predictions in all considered situations agree well with the numerical 
simulations of the complete system. This suggests that the mean field as- 
sumption leading to the kinetic equation should be valid in the considered 
sedimentation model, and the origin of this could be addressed in the future 
with techniques of field theory and renormalization. 

Finally, we have only considered very simple models either without the 
collision efficiency factor, or with a simple forced locality factor conform 
Eq. f fTTj) . Other forms of localizing kernels should be considered for more 
realistic situations. 
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A Locality of power-law distributions 

Power law distributions of the form n{a) ~ a" are important because they 
arise from the formal analysis of the KZ spectra, self-similar solutions, etc. 
However, some of such formal considerations implicitly use convergence of the 
collision integral on RHS of Eq. ([5]) which has the meaning of the interaction 
locality. Conversely, other derivations may assume non-locality i.e. that the 
evolution is dominated mostly by the interactions with the smallest or the 
largest particles in the system corresponding to the vicinities of the small-a 
and the large-cr integration limits. Therefore, the conditions of convergence 
of the collision integral must be found, and this will be done in this appendix 
for a general distribution n{a) ~ a'^. 

Introduce /(cri,(72) = K{ai,a2)nin2. Equation (0) may be expressed in 
terms of / as 

'.(t/2 



dt 



-n 



do-i/((7i,cr - CTi) 



do-i/(o-i,cr) 



We can then split dn/dt into two parts, which we shall call lower and upper 
contributions: 



dt 



-n 



d_ 
dt 



n + 



dt 



n 



with 



d_ 
dt 
d_ 
dt 



(t/2 



n 



n = — 



dcxi[/(cri,cr - ai] 

''max 

dai f{ai,cr) . 



0-/2 



We start by analyzing the lower contribution, more specifically its con- 
vergence as (Jmin goes to 0. For this the value of the integrand at ai ^ a 
needs to be known. This can be approximated by the Taylor expansion 

/((Ti, a - (Ti) - /((Ti, a) ~ aid^f{ai, a) . 

For small cxi we also have f{ai, a) ~ cK^na^^^Uin so we have 



d_ 

dt 



n 



—CK 71 



a/2 



4/3 



n] 
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v<-l 


-l<v<-2 


-2<v 


upper 


local 


non-local 


lower 


non-local 


local 



Table 1: Locality of interaction with small and large particles, as dependent 
on the scaling exponent of n{a) (compare Connaughton et al. [8]). 



The interaction is local at small scales iff the integral above remains finite 
when (Tmin 0. This is equivalent to z/ > — 2. 

We now turn to the upper contribution, more specifically its convergence 
as (Tmax gocs to infinity. For this the value of the integrand at ai ^ cr needs 
to be known. In these asymptotics we have /(ai, a) ~ af^^riin and therefore 



d 
dt 



n ^ —n I nicr^^^do"! 

'a/2 



The interaction is local at large scales iff the integral above remains finite 
when (Tmax oo. This is equivalent to u < —7/3. 

We thus get the picture that for u < —7/3 the interaction is local at 
large scales but non-local at small scales. For —7/3 < z/ < —2 both ends are 
non-local. And for u > —2 interaction is non-local at large scales but local 
at small scales. In particular, we never have locality at both ends. 



B Considerations on the flux 

The flux $(cr) of volume going into particles of volume larger than a can be 
obtained by the following consideration. The flux in question is the volume 
contained in particles of volumes smaller a that merge during unit time with 
some particle to give a particle of volume larger than a. Say one such particle 
has o"i < cr, then it can merge with any particle with a2 such that cri + o"2 > a, 
i.e. a2 > (T — (Ji. Using the collision kernel K the above consideration is made 
formal as 

$(cr) = / dcTi / da2(TiK{ai,a2)n{ai)n{a2) . (29) 

Jo J cr— (71 

One readily verifies by direct computation (and a minor trick) that the right 
hand side of the kinetic equation equals —a~^d„^{cr), so we have as we 
may expect 

. . (30) 
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We immediately remark two things about $. First, it is convergent at the 
lower bound (o"i —>■ 0) if and only if interaction with the small a tail is local, 
and similarly it is convergent at the upper bound (cr2 oo) if and only if 
interaction with the large a tail is local (compare with Appendix 1X1) . 

The other remark is that scales as (j4/3+3+2i/ j^j^ ^ scales as a'^). 

Hence, for u = — (4/3 + 3)/2 = —13/6 we have d„^{cr) = and thus, from 
Eq. (IHUI) . cT~^^/^ is a stationary power law solution. 

The next thing we do is Taylor expand n around ^(cr) in the expres- 
sion (l29l) of the flux. Then to lowest (zeroth) order we get 



with C > (since > 0), and substituting this into Eq. (l30l) we get an 
equation equivalent to Burgers equation (cf. Sect. [7]). 

Remark: perhaps one caveat is that the simple Taylor expansion proposed 
above doesn't seem to correspond to an expansion in some small parameter 
of the problem. One natural small parameter could be g — 1 from the def- 
inition (fTTl) of the collision efficiency. But the expansion in g — 1 would be 
slightly more complex. 

One can carry on this Taylor expansion and get terms of higher order, 
which will have more derivatives d^j and higher powers of a. In the "Burg- 
ers" coordinates introduced in Sect. O the same holds but with powers and 
derivatives in s. In particular to next order we get, in the setup of Sect. 17.21 



dtg = -ds{Cig' + C2sdsg^). 
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